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STATISTICAL  ANALYSIS  OF  THE  OUTPUT  DATA 
FROM  TERMINATING  SIMULATIONS* 

Averill  M.  Law 
University  of  Wisconsin 


ABSTRACT 


In  this  paper  we  precisely  define  the  two  types  of  simulations 
(terminating  and  steady-state)  with  regard  to  analysis  of  simulation 
output  and  discuss  some  common  measures  of  performance  for  each 
type.  In  addition,  we  conclude  from  talking  to  a large  number  of 
simulation  practitioners,  that  a significant  proportion  of  simula- 
tions in  the  real  world  are  of  the  terminating  type.  This  is  con- 
trary to  the  impression  one  gets  from  reading  the  simulation  liter- 
ature, where  the  steady-state  case  is  almost  exclusively  considered. 

Although  analyses  of  terminating  simulations  are  considerably 
easier  than  are  those  of  steady-state  simulations,  they  have  not 
received  a careful  treatment  in  the  literature.  We  discuss  and 
give  empirical  results  for  fixed  sample  size,  relative  width,  and 
absolute  width  procedures  which  can  be  used  for  constructing  con- 
fidence intervals  for  measures  of  performance  in  the  terminating 
case. 


1 

.1 


(j 

; 


* i 


*This  research  was  supported  by  the  Office  of  Naval  Research  under 
contract  N00014-76-C-0403  (NR  047-145) . 


1.  Types  of  Simulations  with  Regard  to  Analysis  of  the  Output 
We  begin  by  giving  a precise  definition  of  the  two  types  of 
simulations  with  regard  to  analysis  of  the  output  data  (cf. 

Gafarian  and  Ancker  [4]  and  Klei jnen  [6]).  A terminating  simulation 
is  one  for  which  any  quantities  of  interest  are  defined  relative 
to  the  interval  of  simulated  time  [O.T^,],  where  T£%  a possibly  de- 
generate random  variable  (r.v.).  Is  the  time  that  a specified  event 
8 occurs.  The  following  are  some  examples  of  terminating  simula- 
tions: 

a)  Consider  a retail  establishment  (e.g.,  a bank)  which  closes  each  even- 
ing (physically  terminating).  If  the  establishment  is  open  from 

9 to  5,  then  the  objective  of  a simulation  might  be  to  estimate 
the  quality  of  customer  service  over  the  8 hour  period.  In  this 
case,  E = (8  hours  of  simulated  time  have  elapsed). 

b)  Consider  a telephone  exchange  which  is  always  open  (physically 
nonterminating).  Since  the  arrival  rate  of  calls  changes  with 
the  time  of  day,  day  of  the  week,  etc.,  it  is  unlikely  that  a 
steady-state  measure  of  performance  (see  Section  2),  which  is 
defined  as  a limit  as  time  goes  to  infinity,  will  exist.  In 
this  case  the  objective  of  a simulation  might  be  to  study  the 
system  during  the  period  of  peak  loading,  say,  of  length  t 
hours,  and  E - {t  hours  of  simulated  time  have  elapsed). 

c)  A system  consists  of  mechanical  and  electronic  components. 


each  of  which  is  subject  to  failure.  The  system  itself  fails 
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when  certain  specified  subsets  of  the  components  fail  and  the 
objective  of  a simulation  might  be  to  estimate  some  character- 
istic of  the  time  to  failure  of  the  whole  system.  In  this  case, 
E = (the  system  fails). 

A steady-state  simulation  is  one  for  which  the  quantity  of 
interest  is  defined  as  a limit  as  the  length  of  the  simulation  goes 
to  infinity.  Since  there  is  no  natural  event  to  terminate  the  simu- 
lation, the  length  of  the  simulation  is  made  large  enough  to  get  a 
"good"  estimate  of  the  quantity  of  interest.  Alternatively,  the 
length  of  the  simulation  could  be  determined  by  cost  considerations. 
The  following  are  some  examples  of  steady-state  simulations: 

a)  Consider  a computer  manufacturer  who  constructs  a simulation 
model  of  a proposed  computer  system.  Rather  than  use  data 

from  the  arrival  process  of  an  existing  computer  system  as  input 
to  the  model,  he  typically  assumes  that  jobs  arrive  in  accord- 
ance with  a Poisson  process  with  rate  equal  to  the  predicted 
arrival  rate  of  jobs  during  the  period  of  peak  loading.  He 
is  interested  in  estimating  the  response  time  of  a job  after 
the  system  has  been  running  long  enough  so  that  initial  condi- 
tions (e.g.,  the  number  of  jobs  in  the  system  at  time  0)  no 
longer  have  any  effect.  (Assuming  that  the  arrival  rate  is 
constant  over  time  allows  steady-state  measures  to  exist.) 

b)  A chemical  manufacturer  constructs  a simulation  model  of  a pro- 
posed chemical  process  operation.  The  process,  when  in 


f 


operation,  will  be  subject  to  randomly  occurring  breakdowns.  The 
input  rate  of  raw  materials  to  the  process  and  the  controllable 
parameters  of  the  process  are  both  assumed  to  be  stationary 
with  respect  to  time.  The  company  would  like  to  estimate  the 
production  rate  after  the  process  has  been  running  long  enough 
so  that  initial  conditions  no  longer  have  any  effect. 

The  remainder  of  this  paper  is  organized  as  follows.  In  Sec- 
tion 2 we  discuss  some  common  measures  of  performance  for  both 
types  of  simulations  and  in  Section  3 we  present  our  findings  on 
the  relative  occurrence  of  each  type  in  the  real  world.  A number 
of  procedures  which  can  be  used  to  construct  confidence  intervals 
for  terminating  simulations  are  discussed  in  Section  4 and,  finally, 
in  Section  5 we  sumnarize  our  findings. 

2.  Measures  of  System  Performance 
To  the  best  of  our  knowledge,  nowhere  in  the  simulation  liter- 
ature are  measures  of  performance  for  terminating  simulations  ex- 
plicitly defined.  In  this  section  we  define  and  contrast  several 
common  measures  of  performance  for  terminating  and  steady-state 
simulations  by  means  of  examples.  (Because  of  the  diversity  of 
terminating  simulations,  it  is  not  possible  to  give  one  definition 
that  fits  all  cases.)  For  the  examples  that  we  consider,  it  is 
possible  to  compute  analytically  measures  of  performance.  Further- 
more, the  same  examples  will  be  considered  in  Section  4 where  we 
discuss  stopping  rules  for  terminating  simulations. 
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A.  Averages 

Consider  first  the  stochastic  process  {Z?;,tal}  for  the  M/M/1 
queue  with  p < 1,  where  D.  is  the  delay  in  queue  of  the  tth  cus- 
tomer.  The  objective  of  a terminating  simulation  of  the  M/M/1 
queue  might  be  to  estimate  the  expected  average  delay  of  the  first 
m customers  given  that  the  number  of  customers  in  the  system  at 
time  0,  .V(0),  is  zero.  The  desired  quantity,  which  we  denote  by 
<i(m|,v(0)=0) , is  then  given  by 


ti(m|/v(0)=0)  = E 


m 

= L £’[^-|A'(0)=0]/m  . 

t= 1 1 


£ D./m |//(0)=0 
i=l  • 


(Alternatively,  if  one  is  interested  in  estimating  the  expected 
average  delay  of  all  customers  who  arrive  and  are  served  in  the 
time  interval  [0,t],  then  the  desired  quantity  is  given  by 


d{t  |.V(0)=0) 


Di/M(t)\N(0)=0 


where  A/(c)  (a  r.v.)  is  the  number  of  customers  who  arrive  and  are 
served  in  the  interval  [0,t].  Note  that  in  this  case  the  expectation 
and  summation  are  not  interchangeable.  Thus,  the  label  "expected 
average  delay"  is  more  general  than  "average  expected  delay".) 

Note  also  that  the  quantity  d(n|jlf(0)s0),  which  is  often  called  a 
transient  characteristic  of  the  stochastic  process  {D.  ,£*1},  explic- 
itly depends  on  the  state  of  the  system  at  time  0;  i.e.,  i(m!Ar(0)=t) 

/ <i(m|iV(0)=j)  for  i 1 j. 


k. 
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The  objective  of  a steady-state  simulation  of  {D^ ,ial}  for 
the  M/M/1  queue  would  be  to  estimate  the  steady-state  expected  av- 
erage delay  d,  which  is  given  by 

d = lim  ii(rrj| jv(0)*t)  for  any  t=0,l,... 
m-*°° 


Observe  that  d is  independent  of  .V(0).  in  Figure  1 we  plot 
J(m|y(0)=0)  as  a function  of  m.  (The  arrival  rate  X = 1 and  the 
service  rate  u = 10/9,  so  p = .9.)  The  horizontal  line  that 
i(rn|y(0)=0)  asymptotically  approaches  is  at  height  d. 

As  a second  example  consider  the  stochastic  process  {£\. ,til } 
for  an  (a,5)  inventory  system  with  zero  delivery  lag  and  backlogging 
where  E.  is  the  expenditure  in  the  ith  period.  This  system  is 
described  in  detail  in  Law  [8].  A possible  objective  of  a terminat- 
ing simulation  would  be  to  estimate  the  expected  average  cost  for 
the  first  m periods  given  that  the  inventory  level  at  the  beginning 
of  period  1 , I, , is  S: 


e(m| Ji=S) 


E^/m\l^S 


The  objective  of  a steady-state  simulation  of  {E. ,tal}  would  be  to 
estimate  the  steady-state  expected  average  cost: 


3 lim  for  any  £=0,±1 ,±2,. . . . 

We  plot  as  a function  of  m and  also  c in  Figure  2. 
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Our  third  example  is  quite  different  from  the  first  two.  In 
the  reliability  model  shown  in  Figure  3 it  is  desired  to  send  a 
message  from  A to  this  will  occur  if  component  1 works  and  either 
component  2 or  3 works.  If  T is  the  time  to  failure  of  the  whole 
system  and  J\  the  time  to  failure  of  component  t(t=l,2,3),  then 

T = min[r^  .inaxCr^  T 3)]  . 

We  further  assume  that  the  Z\' s are  independent  r.v.'s  and  each  has  a 
Weibull  distribution  with  shape  parameter  a=.5  and  scale  parameter  S=l. 

(A  distributional  assumption  for  the  r.'s  is  needed  in  Section  4 where  we 

v 

present  simulation  results  for  this  model.)  The  objective  of  a terminating 
simulation  might  be  to  estimate  the  expected  time  to  failure  of  the 
system  given  that  all  components  are  new,  £(r|al1  components  are  new). 

If  we  assume  that  the  system  is  not  repaired  when  it  fails,  then 
a steady-state  simulation  makes  no  sense  for  this  system.  Such 
could  be  the  case  if  this  system  were  part  of  a space  probe. 

B.  Proportions 

The  usual  criterion  for  comparing  two  or  more  systems  is  some 
sort  of  average  behavior.  However,  different  kinds  of  information 
may  be  of  more  value  in  some  situations.  For  example,  a bank  manager 
may  be  concerned  with  estimating  the  proportion  of  customers  who 
experience  a delay  in  excess  of  5 minutes.  Since  proportions  arp 
really  just  a special  case  of  averages,  we  will  illustrate  them  by 
means  of  the  M/M/1  example. 
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In  a terminating  simulation  of  the  M/M/1  queue  the  objective 
might  be  to  estimate,  instead  of  an  expected  average  delay,  the 
expected  proportion  of  the  first  customers  whose  delay  is  Jess 
than  or  equal  to  x (a  specified  number)  given  that  .V(0)=0.  Denote 
the  desired  quantity  by  P(n,.r  |y(0)"0)  and  let 


1 if  P . { .r 

i 

for  t*l  ,2,... 

|o  if  r.  > .r 


Then  r(m,.r |.V(0)*0)  is  given  by 


r(~!,x|.V(0)°0)  = y 


E V.(.r)A|A’(0)=0 


m 


E i'[V,.(.r)  j A’ ( 0 ) = 0 ] /- 

t»l 


E F{ P . I .V ( 0 ) =0 > /-: 
:*1 


For  a steady-state  simulation,  the  objective  would  be  to  estimate 
the  steady-state  expected  proportion  of  customers  whose  delay  is 
less  than  or  equal  to 


F(x)  = lim  r(nJa*|,v(0)s ) for  any  7=0,1 

3.  Relative  Importance  of  Each  Type  of  Simulation 
Reading  the  simulation  literature  leads  one  to  think  that 
steady-state  simulations  are  more  important;  almost  every  paper 
written  on  the  analysis  of  simulation  output  data  deals  with  the 


kir 


I 


J 


steady-state  case.  This  may  be  a carry-over  from  mathematical 
queueing  theory  where  only  a steady-state  analysis  . > generally 
possible.  However,  we  have  discovered  by  talking  to  a large  number 
of  simulation  practitioners  that  a significant  proportion  of  simula- 
tions in  the  real  world  are  actually  of  the  terminating  type.  The 
following  are  some  reasons  why  a steady-state  analysis  may  not  be 
appropriate: 

a)  The  system  under  consideration  is  physically  terminating.  In 
this  case,  letting  the  length  of  a simulation  be  arbitrarily 
large  makes  no  sense. 

b)  The  input  distributions  for  the  system  change  over  time.  In 
this  case,  steady-state  measures  of  performance  will  probably 
not  exist. 

c)  One  is  often  interested  in  studying  the  transient  behavior  of 
a system  even  if  steady-state  measures  of  performance  exist. 

4.  Stopping  Rules  for  Terminating  Simulations 
In  the  following  three  subsections  we  consider  procedures 
which  can  be  used  to  construct  confidence  intervals  (c.i.'s)  for 
measures  of  performance  for  terminating  simulations.  We  will  not 
consider  the  steady-state  case  since  it  has  been  widely  discussed 
in  the  simulation  literature.  For  surveys  of  fixed  sample  si:e  and 
sequential  procedures  which  can  be  used  to  construct  c.i.'s  for 
steady-state  measures  of  performance , see  Law  [9]  and  law  and 
Kelton  [10],  respectively.  The  random  numbers  used  in  the  remainder 
of  this  paper  were  generated  from  the  generator  discussed  in  [8], 
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A.  Fixed  Sample  Size  Procedures 

Suppose  we  make  n independent  replications  of  a terminating 
simulation.  The  independence  among  replications  is  accomplished 
by  using  different  random  numbers  for  each  replication  and  by 
starting  each  one  with  the  same  initial  conditions.  If  X.  is  the 
estimator  of  interest  from  the  ;'th  replication  ( t=l  ,2 , . . . ,n) , then 
the  .V . ' s are  independent  identically  distributed  (i.i.d.)  r.v.'s. 

(For  the  M/M/1  queue,  l . might  be  the  average  7 P ./n  or  the  pro- 

m k i s ] * 

* i 1 

portion  7 Y.(x)/m.)  If,  in  addition,  the  .Y.'s  are  normally  dis- 

j-1  J 1 

tributed,  then  a 100(1  -*)%  (0<a<l)  c.i.  for  u = F(.Y)  is  given  by 

*(»)  * 

where  X(n)  and  ?2  {>.)  are  the  usual  sample  mean  and  variance,  re- 
spectively, and  tf  ^ is  the  1 “ «/2  point  for  a t distribu- 

tion with  •;-!  degrees  of  freedom. 

In  practice  the  .Y.'s  will  not  be  normally  distributed  and  the 
c.i.  given  by  (1)  will  be  only  approximate.  To  investigate  the 
effect  of  nonnormality,  we  simulated  the  ♦'hree  stochastic  models 
of  Section  2.  For  the  M/M/1  queue  with  p = .9,  the  (c  ,S)  inven- 
tory system,  and  the  reliability  model,  respectively,  the  quantities 
of  interest  were  d(2b  |ff(0)=0)  = 2.12,  o(12  | = 99.52,  and 

F(7'|all  components  are  new)  = .778.  (See  [81  fora  discussion 
of  how  to  compute  the  first  two  quantities.)  For  each  model  we 
performed  500  independent  simulation  experiments,  for  each  experi- 
ment we  considered  n = 5,  10,  20,  40.  and  for  each  n we  used  (1) 


- 10  - 


to  construct  a 90%  c.i.  for  the  desired  quantity.  In  Tables  1,2, 
and  3 we  give  the  proportion,  p,  of  the  500  c.i.'s  which  covered  the 
desired  quantity,  a 90%  c.i.  for  the  true  coverage,  and  the  average 
value  of  the  c.i.  half  length  divided  by  the  point  estimate  over 
the  500  experiments  for  the  three  models.  The  90%  c.i.  for  the 
true  coverage  was  computed  from 

p ± 1 .645/p(l-p)/500  . 

Observe  that  for  the  M/M/1  queue  and  the  ( s ,S ) inventory  system  the 
coverages  are  quite  close  to  90%,  but  for  the  reliability  model 
there  is  a significant  degradation  in  coverage,  apparently  caused 
by  a severe  departure  from  normality.  To  see  if  this  is  indeed 
the  case,  we  generated  1000  x. 's  for  each  stochastic  model  and 
estimated  the  skewness  and  kurtosis.  These  estimates,  which  are 
presented  in  Table  4,  indicate  that  the  X.'s  for  the  reliability 
model  are  considerably  more  nonnormal  than  are  those  for  the  other 
two  models.  This  conclusion  was  reinforced  by  plotting  histograms 
for  the  three  sets  of  data. 

B.  Relative  Width  Procedures 

One  disadvantage  of  the  fixed  sample  size  approach  to  construct- 
ing a c.i.  is  that  the  simulator  has  no  control  over  the  c.i.  half 
length;  for  fixed  n,  the  half  length  will  depend  on  the  population 
variance  a2  = Var(x).  In  this  subsection  we  consider  two  sequential 
procedures  which  allow  one  to  specify  the  "relative  precision"  of 


■u,n **• 


- 11  - 


a c.1.  Both  assume  that  A^.A^,...  is  a sequence  of  i.i.d.  r.v.'s 
which  need  not  be  normal. 

The  first  procedure  has  been  suggested  for  use  In  several 
different  contexts;  see  Iglehart  [5],  Lavenberg  and  Sauer  [7],  and 
Thomas  [13].  The  objective  of  the  procedure  is  to  construct  a 
100(1 -ot)%  c.i.  for  u such  that  the  difference  between  the  point 
estimator  X(n)  and  p is  no  more  than  100  y%  of  X(n) , that  is, 

|x(n)-p|  if  y | A'(n)  | for  0<y<l  • (2) 


Choose  an  initial  sample  size  * 2,  let 


and  let 


Nr,  l(Y,a)  = min|n:  n * nQ>  « («)  > °»  -j"  ^ Y>  • 


(3) 


(Note  that  /V(>  ^y.a),  which  is  the  required  number  of  replications, 
is  a r. v. ) Then  use 


7r,l( Y,ot)  = 1 (y  .a)  (iVr  j (y.«)  •«)] 


(4) 


as  an  approximate  100  (l-u)T.  c.i.  for  p.  It  easily  follows  from 
(3)  and  (4)  that  1 ,(y,oi)  satisfies  the  criterion  given  by  (2). 

X*  9 I 

Furthermore,  using  an  argument  similar  to  the  one  employed  by 
Lavenberg  and  Sauer  in  the  context  of  the  regenerative  method  for 
steady-state  simulations,  we  have  been  able  to  prove  the  following 


theorem. 


r 
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Theorem  1.  If  u f 0 and  0 < az  < then  11m  P{pC  T , (y.<*)) 

rK)+  r**1 

3 1 - a. 

The  objective  of  the  second  procedure,  which  is  due  to  Nadas 
[11],  is  to  construct  a c.i.  such  that 

|.Y(«)-vi|  *•  y |u|  for  0 < y < 1 . (5) 

Let 

I ’2{n)  3 jl  + [.VJ.-Y(n)];'l/«  = (1/m)  + (h-1)s’(m)/m  . 

5n,2<"’a)  * tn.},}-a/2^vlM/n  ’ 

and 

\ <5„  2(”*a)  ) 

Nr,2^’a)  * ■1Y:»  * V "feT  v<  y\  - 

Then  use 

^,2(Y.a)  « 

as  an  approximate  lOO(l-a);.  c.i.  for  u.  From  (6)  it  is  easy  to 
show  that  ,,(>,a)  satisfies  the  criterion  given  by  (5).  Further- 
more, the  following  theorem  was  proved  by  Nadas. 

Theorem  2.  If  u / 0 and  0 < o'1  < «\  then  lim  /’(u  C I ?(y»^)1 

r*o  r* 

= 1 - a. 

In  order  to  compare  the  two  procedures  and  to  determine  the 
effect  of  non-infinitesimal  y on  coverage,  we  once  again  simulated 
the  three  stochastic  models.  For  each  model  we  performed  500  inde- 
pendent experiments,  for  the  M/M/1  queue  and  the  reliability  model 


v(.Y  , 
r,2 


, •)  ( V * ll* ) ) 

fvy- 


.(y,u)  ) 


I 


13  - 


we  considered  y = .2,  .1,  .05  for  each  experiment,  and  for  the  In- 
ventory system  we  considered  y ■ .2,  .1,  .05,  .025,  .0125,  .00625 
for  each  experiment.  In  all  cases,  rig  » 5.  In  Tables  5,  6,  and  7 
we  give  point  estimates  and  90*  c.i.'s  for  the  true  coverages,  point 
estimates  and  90*  c.i.'s  for  K{N  t . (y.u)  }(i'=l  ,2) , and  the  average 
c.1.  half  lengths  over  the  500  experiments.  We  considered  more 
values  of  y for  the  inventory  system  because  It  appeared  from  our 
empirical  results  that  a smaller  y is  required  for  the  coverage  ul- 
timately to  converge  to  the  desired  level.  (A  smaller  y is  required 
for  this  model  to  get  a large  value  of  N .(y,a).)  Note  also  that 

V ft* 

convergence  of  coverage  does  not  appear  to  be  monotone. 

We  repeated  the  above  500  experiments  using  the  same  random 
numbers  and  Nq  = 2.  For  procedure  2 the  results  were  identical; 
however,  for  procedure  1 there  was  a significant  degradation  in 
coverage  due  to  premature  stopping  on  replications  2,3,  or  4.  For 
example,  the  coverage  for  the  M/M/1  queue  with  y = .2  was  .798. 

c.  Absolute  Width  Procedures 

In  this  subsection  we  present  two  procedures  which  allow  one 
to  construct  a 100(l-a)*  c.1.  for  u such  that 

|*(n>-u|  < * , (7) 

where  <•  is  a specified  positive  number. 

The  first  procedure,  which  is  due  to  Chow  and  Robbins  [1],  as- 
sumes that  is  a sequence  of  i.i.d.  r.v.'s.  Choose  * 2. 
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Let  v2(n)  be  defined  as  in  Subsection  4.B,  let 


N ,(c,ot)  = minium  a n , y2(n)  £ 77 — y 

a']  0 (Vl,l-a/2' 


and  then  use 


I . [a, a)  = [X(/V  ,{at a))  - a,X(N  ,(c,a))  + <?] 

<2,1  u , 1 u,  1 

as  an  approximate  100(l-ot)%  c.i.  for  y.  It  is  clear  that  Iq  ^(c,a) 
satisfies  the  criterion  given  by  (7).  The  following  theorem  was 
proved  by  Chow  and  Robbins. 

Theorem  3.  If  0 < a2  < °°,  then  lim  P{y  C -T  -.(c.ct)  = 1 - a. 

<r*C+  a’ 

For  an  empirical  evaluation  of  the  above  procedure  under  the  assump 
tion  that  the  X.'s  are  normal,  see  Starr  [12]. 

The  second  procedure,  which  is  due  to  Dudewicz  [2],  assumes 
that  the  X/s  are  i.i.d.  normal  r.v.'s.  Initially  make  »qUq*2) 
replications  of  the  simulation  and  compute  x(rcg)  and  s2(«g).  Let 

Na  2^a,cx^  = Hiax{n0+l,rW2S2(n0)l}  » 


where  u = t , , /0/c*  and  fsl  is  the  smallest  integer  ^ z.  Make 

«g-l,l-a/2 

N 0(c,a)  - additional  replications  of  the  simulation,  let 
a,  i U 

2 (c7',ot) 

Hn  2(^,a)-^0)  = £ Xi^Na  2^°^~n0^* 

* i=n0+1 

and  let  X{Nq  2(c?,a))  = a^ng)  + a2?(lVa  2(<?,a)-ng),  where 


al  ‘ N 2(c?,a) 
a%c 


1 + 


V ■ »0  V “Js  ("o1  L 
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and  u2  = 1 - u1 . Then  use 

I = [A"U  Ac,  <x))-c,X[N  ?(o,ot))+c] 

a,t  cx,c  u, c 

as  an  approximate  100(l-a)%  c.i.  for  p.  Dudewicz  has  proved  the 
following  theorem. 

Theorem  4.  P{\i  C I ?(t>,a)}  = 1 - a for  all  a > 0. 

To  compare  the  sequential  procedure  of  Chow  and  Robbins  and 
the  two-stage  procedure  of  Dudewicz,  we  performed  500  independent 
experiments  for  each  model.  To  make  the  absolute  width  results 
somewhat  comparable  to  the  relative  width  results,  we  chose  the 
values  of  to  correspond  to  the  values  of  y;  that  is,  for  each  y 
we  chose  a - yp.  For  the  Chow  and  Robbins  procedure  we  chose 

= 5 and  for  the  Dudewicz  procedure  we  considered  = 15,  30,  and 
60.  (Dudewicz  [3]  recommended  that  nQ  be  at  least  12.)  The  re- 
sults of  the  simulation  experiments  for  the  three  models  are  given 
in  Tables  8,  9,and  10. 

5.  Summary  and  Conclusions 

We  have  defined  terminating  and  steady-state  simulations  and 
have  discussed  some  common  measures  of  performance  for  each  type. 

In  addition,  we  have  concluded  from  talking  with  simulation  practi- 
tioners that  a significant  proportion  of  real-world  simulations  are 
of  the  terminating  type.  This  is  fortunate  because  it  means  that 
classical  statistical  analysis  for  i.i.d.  observations  (e.g.,  con- 
fidence intervals,  hypothesis  testing,  ranking  and  selection,  etc.) 


I 


I 
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is  applicable  to  analyzing  many  simulations.  On  the  other  hand, 
in  the  steady-state  case  there  is  still  not  a totally  acceptable 
procedure  even  for  the  relatively  simple  problem  of  constructing 
a c.i.  for  a steady-state  expected  average. 

We  have  also  considered  procedures  for  constructing  c.i.'s 
for  terminating  simulations.  If  one  is  performing  an  exploratory 
experiment  where  precision  of  the  c.i.  may  not  be  overwhelmingly 
important,  then  we  recommend  using  a fixed  sample  size  procedure. 
However,  if  the  A\'s  are  highly  nonnormal  and  if  the  number  of  rep- 
lications n is  too  small,  then  the  actual  coverage  of  the  con- 
structed c.i.  may  be  considerably  lower  than  that  desired  (see 
Table  3). 

If  one  wants  a c.i.  having  half  length  that  is  small  relative 
to  the  point  estimate,  then  a relative  width  procedure  may  be  used. 
We  recommend  using  Procedure  2 (due  to  Nadas)  With  «g  * 5.  Proce- 
dure 2 appears  to  give  slightly  better  coverage,  its  criterion  (see 
(5))  is  more  intuitive  than  the  criterion  of  Procedure  1 (see  (2)), 
and  Procedure  2 does  not  seem  subject  to  premature  stopping  even 
for  Nq  = 2.  (On  the  other  hand.  Procedure  1 uses  a more  intuitive 
expression  to  construct  a c.i.) 

If  one  wants  a c.i.  for  which  the  half  length  is  a specified 
number,  then  an  absolute  width  procedure  may  be  used.  We  recommend 
using  the  Chow  and  Robbins  procedure  with  Mg  * 5.  Their  procedure 
generally  requires  a smaller  average  sample  size,  the  variance  of 
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the  sample  size  is  smaller,  and  its  coverage  seems  to  be  less  af- 
fected by  departures  from  normality  (see  Table  10). 

In  general,  we  believe  that  relative  width  procedures  are  more 
useful  than  absolute  width  procedures  due  to  the  difficulty  in 
specifying  the  absolute  width  a for  most  simulation  experiments. 

When  using  either  the  Nadas  procedure  or  the  Chow  and  Robbins  pro- 
cedure, we  believe  that  it  is  advisable  to  choose  a y or  a which 
will  cause  the  procedure  to  run  until  the  sample  is  at  least  of  mod- 
erate size;  perhaps,  at  least  30.  (Since  both  procedures  are  based 
on  the  central  limit  theorem,  it  is  unreasonable  to  think  that  they 
will  work  well  in  general  for  a small  sample  size;  see  the  results 
for  y = .025  in  Table  6.)  Finally,  we  mention  that  precise  c.i.'s 
may  be  unaffordable  in  the  real  world  due  to  the  high  cost  of  mak- 
ing a single  replication. 
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Table  4.  Skewness  and  kurtosis  for  the  Three 

Stochastic  Models  and  the  Normal  Distribution. 


Stochastic  Model 
or  Distribution 

Skewness 

Kurtosis 

Normal  Distribution 

0* 

3 * 

: M/M/1  Queue 

1 .66 

6.43 

! (s,S)  Inventory  System 

.45 

3.76 

j Rel labi 1 1 ty  Model 

5.18 

54.39 

♦Theoretical  Values 


Table  5.  Relative  Width  Results  for  j(25 |A'(0)=0)  = 2.12,  M/M/1 
Queue  with  p = .9. 


Procedure  1 

Procedure 

D 

b'{Np  j (y, a)}  j coverage 

— 

average  c. i . 
half  length 

E { N ^ £ ( Y ) 

coverage 

average  c . i . 
half  length 

— 

.2 

42.3+0.9 

i .842 ±.02 7 

.414 

.862+. 025 

.437 

.1 

175.2+1 .7 

.860±.026 

.211 

174.5±1.7 

.868+. 025 

.213 

.05 

704.4+3.5 

j ,884± .024 

703.7+3.5 

.882 +.024 

.106 

Table  6.  Relative  Width  Results  for  ^(121^=5)  = 99.52,  (s,S)  Inventory 
System. 


Procedure  1 

r 

Procedure  2 

fl 

RE 

coverage 

D 

coverage 

average  c.i. 
half  length 

.2 

.902+. 022 

4.89 

5. 0+0.0 

20.74 

.1 

5. 0+0.0 

.902+. 022 

4.86 

5.0±0.0 

.05 

5. 9+0.1 

.892+. 023 

3.97 

5 . 7±0 . 1 

.962+. 014 

4.99 

.025 

13.3+0.4 

.834+. 027 

2.35 

12.3+0.4 

.858+.C2C 

2.48 

.0125 

51.0+1.0 

.856+. 026 

1.23 

49.8+1.0 

1.24 

.00625 

206.3+1.8 

.872+. 025 

0.62 

205. 4±1 .8 

0.62 



Table  7.  Relative  Width  Results  for  £"(!T | al  1 components  new)  = .778, 
Rel iabi 1 ity  Model . 


Procedure  1 

Procedure  2 

I 

B{Np^(yta)} 

coverage 

average  c.i. 
half  length 

k{^2(y,a)} 

! | 

coverage  j 

21 3.7±4.5 

.876+. 024 

BH 

214.1+4.5 

.908+. 021 

ihB 

907.4+11 .2 

.898 +.02 2 

E3 

908.6+10.8 

.902+. 022 

3720.5+23.7 

.882+. 024 

.039 

3720.0+23.7 

.884±.024 

Table  8.  Absolute  Width  Results  for  d(25 |tf(0)=0)  = 2.12,  M/M/1 
Queue  with  p = .9. 


Chow  and 

Robbins 

Dudewicz 

o 

iTv  TCT-TOT 

U j 1 

coverage 

7!0 

TO,  Aota)} 

'coverage 

.425 

L ■ — -T- 

38. Oil. 2 

.800±.029 

15 

30 

60 

49.9i2.1 
48.2il .3 

62 . 1*0. 4 

.850*  .026 
.912i .020 
.926*  .019 

.212 

173.5±2.5 

.89&t,022 

1 5 

30 

60 

m>8i 

185.7*5.6 

182.9*4.0 

. 854*  .026' 
.888*  .023 
.8941.023 

.106 

706.8±4.8 

. 906i . 02 1 

15' 

30 

60 

786.1*34.2 

741.1*22.6 

730.2*15.7 

.868* .025 
.878-1.024 
. 8981 . 022 

Table  9.  Absolute  Width  Results  for  ^(121/,=^)  = 99.52,  U*,5) 
Inventory  System. 


" “Chow  anJToWins'T  D'udewl cY  ~ " I 

r~ 

E(N^  ^ (c*, a) ) (coverage  j>;0 

y{X 

-UJ  c _ 

coverage 

119.90 

5.010.0 

1 1 5 

1.0  30 

6° 

16.010.0 

31.010.0 

61 .0*0.0 

.936* .018 
.8781.024 
.8881.023 

r 

9.95 

5. 0*0.0 

J 

15H 

1.0  j 30 

1 60 

16.010.0 

31.01-0.0 

61.0*0.0 

.936* .018 
.8801.024 
.8901.023 

4.98 

5.710.1 



.9761.011 

30 

[60 

rT6:T)  *0.0  ■ H 
31.0*0.0 
61.0*0.0 

.8821.024 

.8861.023 

2.49 

" 

12.310.4 



.8801.024 

IT 

30 

601 

18.5*0.3 

31.010.0 

61.0*-0.0 

79081.021 

.8941.023 

.8821.024 

1.24 

48. 3H.  1 



' 

.8721.025 

'15 

30 

60 

60.8*270 

55. OH. 3 
62.910.4 

.9041.022 

.9121.020 

.9021.022 

0.62 

204.4*1.8 

.8961.022 

15" 

30 

60 

241 .718.1 

217.915.1 
211.513.4 

.9121.020 

.8981.022 

.9121.020 

c 

.156 

.078 

.039 


Chow  and 
E{  N Ac,  a) 

'±JL1 

1 79 . 5±7 . 0 


888. 0± 14. 5 


3672.1 ±32. 9 


Robbins 


coverage 


. 774± .031 


.900+.022 


.884± .024 


Dudewicz 


coverage 


.7 

220. 8±  17.0  .7 
231. 7±  14.9  .8 


880. 6±  68.0  .794±.030 
922. 2±  59.6  .838+.027 


3925.6±435.8 

3520.91272.0 

3687.21238.5 


Figure  2.  e(m|l-|=a)  as  a Function  of  m for  the  (e,S)  Inventory 

System. 
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